Exercise 3: electron vs gamma separation¶
1. Introduction¶
Fig. 14 Difference between an electron vs photon is at the start of the electromagnetic shower, where the photon has a gap. From Wouter Van De Pontseele, ICHEP 2020.¶
Electrons are visible in a LArTPC detector because of the electromagnetic showers that they trigger.
Photons, on the other hand, are neutral (no charge) and thus remain invisible to the LArTPC eyes until they convert into electrons (pair production) or Compton scatter. In both cases, the visible outcome will be an electromagnetic shower.
How can we differentiate the two, then? The answer is in the very beginning of the EM shower. For an electron, this shower will be topologically connected to the interaction vertex where the electron was produced. For a photon, there will be a gap (equal to the photon travel path) until the EM shower start (when the photon becomes indirectly visible through pair production or Compton scatter). That seems simple enough, right? Wrong, of course.
Energetic photons could interact at a distance short enough from the interaction vertex, that we would not be able to see the gap. Or, the hadronic activity might be invisible, because it includes neutral particles or because the particles are too low energy to be seen. In that case the interaction vertex might be hard to identify, and the notion of a gap goes away too. For such cases, fortunately, there is another way to tell electrons from gamma showers. Another major difference is in the energy loss rate at the start of the EM shower. An electron would leave ionization corresponding to a single ionizing particle, whereas a pair of electron + positron coming from a photon pair production would add up to two ionizing particle. Thus, we expect the dE/dx at the beginning of the shower to be roughly twice larger in the case of a gamma-induced shower compared to an electron-induced shower.
Fig. 15 Example from MicroBooNE. Left is the shower \(dE/dx\), right is the gap between the vertex and shower start. From Wouter Van De Pontseele, ICHEP 2020.¶
Why do we care? The difference becomes significant if, for example, you are looking for electron neutrinos. One of the key signatures you would be looking for are electrons.
In this exercise, we will focus on finding the start of EM showers and computing the reconstructed dQ/dx in these segments. Optionally, you could compare that to the result of using automatic PID as predicted by the chain.
2. Setup¶
a. Software and data directory¶
import os, sys
SOFTWARE_DIR = '%s/lartpc_mlreco3d' % os.environ.get('HOME')
DATA_DIR = os.environ.get('DATA_DIR')
# Set software directory
sys.path.append(SOFTWARE_DIR)
b. Numpy, Matplotlib, and Plotly for Visualization and data handling.¶
import numpy as np
import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
'figure.figsize':(15, 10),
})
seaborn.set_context('talk')
import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)
c. MLRECO specific imports for model loading and configuration setup¶
from mlreco.main_funcs import process_config, prepare
import warnings, yaml
warnings.filterwarnings('ignore')
cfg = yaml.load(open('%s/inference.cfg' % DATA_DIR, 'r').read().replace('DATA_DIR', DATA_DIR),Loader=yaml.Loader)
process_config(cfg, verbose=False)
/usr/local/lib/python3.8/dist-packages/MinkowskiEngine/__init__.py:36: UserWarning:
The environment variable `OMP_NUM_THREADS` not set. MinkowskiEngine will automatically set `OMP_NUM_THREADS=16`. If you want to set `OMP_NUM_THREADS` manually, please export it on the command line before running a python script. e.g. `export OMP_NUM_THREADS=12; python your_program.py`. It is recommended to set it below 24.
Config processed at: Linux tur024 3.10.0-1160.42.2.el7.x86_64 #1 SMP Tue Sep 7 14:49:57 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
$CUDA_VISIBLE_DEVICES="0"
d. Initialize and load weights to model using Trainer.¶
# prepare function configures necessary "handlers"
hs = prepare(cfg)
dataset = hs.data_io_iter
Welcome to JupyROOT 6.22/09
Loading file: /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/mpvmpr_012022_test_small.root
Loading tree sparse3d_reco
Loading tree sparse3d_reco_chi2
Loading tree sparse3d_pcluster_semantics_ghost
Loading tree cluster3d_pcluster
Loading tree particle_pcluster
Loading tree particle_mpv
Loading tree sparse3d_pcluster_semantics
Loading tree sparse3d_pcluster
Loading tree particle_corrected
Warning in <TClass::Init>: no dictionary for class larcv::EventNeutrino is available
Warning in <TClass::Init>: no dictionary for class larcv::NeutrinoSet is available
Warning in <TClass::Init>: no dictionary for class larcv::Neutrino is available
Shower GNN: True
Track GNN: True
Particle GNN: False
Interaction GNN: True
Kinematics GNN: False
Cosmic GNN: False
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Shower GNN: True
Track GNN: True
Particle GNN: False
Interaction GNN: True
Kinematics GNN: False
Cosmic GNN: False
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Ghost Masking is enabled for UResNet Segmentation
Restoring weights for from /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/weights_full_mpvmpr_012022.ckpt...
Done.
Let’s load one iteration worth of data into our notebook:
data, result = hs.trainer.forward(dataset)
Segmentation Accuracy: 0.9824
PPN Accuracy: 0.8162
Clustering Accuracy: 0.9279
Shower fragment clustering accuracy: 0.9255
Shower primary prediction accuracy: 0.7553
Track fragment clustering accuracy: 0.9480
Interaction grouping accuracy: 0.9705
Particle ID accuracy: 0.8167
Primary particle score accuracy: 0.8810
e. Setup Evaluator¶
from analysis.classes.ui import FullChainEvaluator
# Only run this cell once!
evaluator = FullChainEvaluator(data, result, cfg, deghosting=True)
print(evaluator)
FullChainEvaluator(num_images=10)
entry = 4 # Batch ID for current sample
print("Batch ID = ", evaluator.index[entry])
Batch ID = 4
3. Identifying Shower Primaries¶
Step 1: Get shower primary fragments¶
By using the only_primaries=True option, we can select out primary particles in this image. We will also load true_particles for comparison.
particles = evaluator.get_particles(entry, only_primaries=True)
true_particles = evaluator.get_true_particles(entry, only_primaries=True)
from pprint import pprint
pprint(particles)
[Particle( Image ID=0 | Particle ID=0 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.09% | Interaction ID: 3 | Size: 1026 ),
Particle( Image ID=0 | Particle ID=2 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.52% | Interaction ID: 3 | Size: 63 ),
Particle( Image ID=0 | Particle ID=3 | Semantic_type: Shower Fragment | PID: Electron | Primary: 1 | Score = 97.70% | Interaction ID: 3 | Size: 1509 )]
Alternatively, as you may have noticed, the primariness information is also stored in the Particle instance as an attribute with name is_primary. If you prefer to view the full image and then select out primaries manually:
particles = evaluator.get_particles(entry, only_primaries=False)
true_particles = evaluator.get_true_particles(entry, only_primaries=False)
from pprint import pprint
pprint(particles)
[Particle( Image ID=0 | Particle ID=0 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.09% | Interaction ID: 3 | Size: 1026 ),
Particle( Image ID=0 | Particle ID=1 | Semantic_type: Shower Fragment | PID: Photon | Primary: 0 | Score = 56.53% | Interaction ID: 1 | Size: 12 ),
Particle( Image ID=0 | Particle ID=2 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.52% | Interaction ID: 3 | Size: 63 ),
Particle( Image ID=0 | Particle ID=3 | Semantic_type: Shower Fragment | PID: Electron | Primary: 1 | Score = 97.70% | Interaction ID: 3 | Size: 1509 ),
Particle( Image ID=0 | Particle ID=4 | Semantic_type: Track | PID: Muon | Primary: 0 | Score = 67.55% | Interaction ID: 4 | Size: 1926 ),
Particle( Image ID=0 | Particle ID=5 | Semantic_type: Track | PID: Pion | Primary: 0 | Score = 34.94% | Interaction ID: 1 | Size: 197 ),
Particle( Image ID=0 | Particle ID=6 | Semantic_type: Delta Ray | PID: Electron | Primary: 0 | Score = 34.86% | Interaction ID: 4 | Size: 23 ),
Particle( Image ID=0 | Particle ID=7 | Semantic_type: Delta Ray | PID: Electron | Primary: 0 | Score = 84.00% | Interaction ID: 1 | Size: 25 ),
Particle( Image ID=0 | Particle ID=8 | Semantic_type: Delta Ray | PID: Electron | Primary: 0 | Score = 31.53% | Interaction ID: 4 | Size: 21 )]
Let’s quickly plot the particles and visualize which ones are predicted as primaries. Here is one way to do it with the trace_particles function:
from mlreco.visualization.plotly_layouts import white_layout, trace_particles, trace_interactions
traces = trace_particles(particles, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_true = trace_particles(true_particles, color='is_primary', colorscale='rdylgn')
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
The green voxels are predicted primary particles, while red indicates non-primary.
It is often easier to further break down the shower into different fragments and locate which of the shower fragments actually correspond to a predicted primary.
fragments = evaluator.get_fragments(entry)
traces = trace_particles(fragments, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_right = trace_particles(fragments, color='id', colorscale='rainbow') # This time, we'll plot the predicted particle
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_right, rows=[1] * len(traces_right), cols=[2] * len(traces_right))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
# TODO: Plot true fragment labels
Step 2: Identify the startpoint of the shower primary¶
During initialization of the Particle instance, PPN predictions are assigned to each particle if the distance between then is less than a predetermined threshold (attaching_threshold). PPN predictions that are matched to particles in this way are then stored in each Particle instance as attributes (ppn_candidates)
print("Minimum voxel distance required to assign ppn prediction to particle fragment = ", evaluator.attaching_threshold)
Minimum voxel distance required to assign ppn prediction to particle fragment = 2
fragments = evaluator.get_fragments(entry, only_primaries=False)
The first three columns are the \((x,y,z)\) coordinates of the PPN points. The fourth column is the PPN prediction score, and the last column indicates the predicted semantic type of the point.
We first visualize whether the predicted ppn candidates accurately locate the shower fragment start:
traces = trace_particles(fragments, color='id', size=1, scatter_ppn=True, highlight_primaries=True) # Set scatter_ppn=True for plotting PPN information
traces_true = trace_particles(true_particles, color='id', size=1)
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
The left scatterplot highlighits primary shower fragments and its ppn candidates, while non-primaries are showed with faded color. The right plot shows true particle labels.
Identifying the primary shower fragments (as above) allow us to select all the voxels of the primary fragment which are close to the shower start, i.e. within some radius of the predicted PPN shower point. Of course, as expected from the scatterplot above, we may also include some cuts on the total voxel count to pick shower primary fragments that are large enough for our \(dQ/dx\) analysis.
For convenience, from now on we will only work with primary fragments:
fragments = evaluator.get_fragments(entry, only_primaries=True)
Step 3. Compute \(dQ/dx\) near the shower start¶
Let’s first fix some parameters for our \(dQ/dx\) computation. Let’s say the we select all points within a radius of 10 voxels from the predicted PPN shower start point of a given primary fragment and require that the selected segment size should at least be 3 voxels long.
from sklearn.decomposition import PCA
from scipy.spatial.distance import cdist
min_segment_size = 3 # in voxels
radius = 10 # in voxels
pca = PCA(n_components=2)
Write a compute_shower_dqdx function that takes a list of primary fragments and returns a list of computed dQ/dx values for each fragment.
def compute_shower_dqdx(frags, r=10, min_segment_size=3):
'''
Inputs:
- frags (list of ParticleFragments)
Returns:
- out: list of computed dQ/dx for each fragment
'''
out = []
for frag in frags:
assert frag.is_primary # Make sure restriction to primaries
if (frag.startpoint < 0).any():
continue
ppn_prediction = frag.startpoint
dist = cdist(frag.points, ppn_prediction.reshape(1, -1))
mask = dist.squeeze() < r
selected_points = frag.points[mask]
if selected_points.shape[0] < 2:
continue
proj = pca.fit_transform(selected_points)
dx = proj[:, 0].max() - proj[:, 0].min()
if dx < min_segment_size:
continue
dq = np.sum(frag.depositions[mask])
out.append(dq / dx)
return out
compute_shower_dqdx(fragments)
[566.5370439012097, 8604.684675385452, 3761.2255977052623, 4063.2545505781763]
Step 4. Collect data over multiple images and plot results¶
iterations = 10
collect_dqdx = []
for iteration in range(iterations):
data, result = hs.trainer.forward(dataset)
evaluator = FullChainEvaluator(data, result, cfg, deghosting=True)
for entry, index in enumerate(evaluator.index):
# print("Batch ID: {}, Index: {}".format(entry, index))
fragments = evaluator.get_fragments(entry, only_primaries=True)
dqdx = compute_shower_dqdx(fragments, r=radius, min_segment_size=min_segment_size)
collect_dqdx.extend(dqdx)
collect_dqdx = np.array(collect_dqdx)
Segmentation Accuracy: 0.9851
PPN Accuracy: 0.8156
Clustering Accuracy: 0.9136
Shower fragment clustering accuracy: 0.8993
Shower primary prediction accuracy: 0.6915
Track fragment clustering accuracy: 0.9334
Interaction grouping accuracy: 0.9381
Particle ID accuracy: 0.8030
Primary particle score accuracy: 0.9083
Segmentation Accuracy: 0.9688
PPN Accuracy: 0.7979
Clustering Accuracy: 0.8609
Shower fragment clustering accuracy: 0.9100
Shower primary prediction accuracy: 0.6218
Track fragment clustering accuracy: 0.9569
Interaction grouping accuracy: 0.8885
Particle ID accuracy: 0.8219
Primary particle score accuracy: 0.8852
Segmentation Accuracy: 0.9932
PPN Accuracy: 0.8321
Clustering Accuracy: 0.9033
Shower fragment clustering accuracy: 0.9487
Shower primary prediction accuracy: 0.7103
Track fragment clustering accuracy: 0.9720
Interaction grouping accuracy: 0.9083
Particle ID accuracy: 0.8269
Primary particle score accuracy: 0.9314
Segmentation Accuracy: 0.9733
PPN Accuracy: 0.7615
Clustering Accuracy: 0.8377
Shower fragment clustering accuracy: 0.9361
Shower primary prediction accuracy: 0.5484
Track fragment clustering accuracy: 0.9534
Interaction grouping accuracy: 0.9403
Particle ID accuracy: 0.7667
Primary particle score accuracy: 0.9580
Segmentation Accuracy: 0.9814
PPN Accuracy: 0.7792
Clustering Accuracy: 0.9067
Shower fragment clustering accuracy: 0.9515
Shower primary prediction accuracy: 0.6935
Track fragment clustering accuracy: 0.9500
Interaction grouping accuracy: 0.9265
Particle ID accuracy: 0.8947
Primary particle score accuracy: 0.9134
Segmentation Accuracy: 0.9819
PPN Accuracy: 0.8271
Clustering Accuracy: 0.8859
Shower fragment clustering accuracy: 0.8776
Shower primary prediction accuracy: 0.5574
Track fragment clustering accuracy: 0.9530
Interaction grouping accuracy: 0.9259
Particle ID accuracy: 0.8868
Primary particle score accuracy: 0.9298
Segmentation Accuracy: 0.9840
PPN Accuracy: 0.8182
Clustering Accuracy: 0.8362
Shower fragment clustering accuracy: 0.8622
Shower primary prediction accuracy: 0.7386
Track fragment clustering accuracy: 0.9604
Interaction grouping accuracy: 0.9454
Particle ID accuracy: 0.9091
Primary particle score accuracy: 0.8763
Segmentation Accuracy: 0.9723
PPN Accuracy: 0.8146
Clustering Accuracy: 0.8628
Shower fragment clustering accuracy: 0.7710
Shower primary prediction accuracy: 0.6327
Track fragment clustering accuracy: 0.8976
Interaction grouping accuracy: 0.9328
Particle ID accuracy: 0.7045
Primary particle score accuracy: 0.7976
Segmentation Accuracy: 0.9494
PPN Accuracy: 0.8093
Clustering Accuracy: 0.8473
Shower fragment clustering accuracy: 0.8637
Shower primary prediction accuracy: 0.4839
Track fragment clustering accuracy: 0.9450
Interaction grouping accuracy: 0.9308
Particle ID accuracy: 0.8000
Primary particle score accuracy: 0.8991
Segmentation Accuracy: 1.0000
PPN Accuracy: 0.6214
Clustering Accuracy: 0.9372
Shower fragment clustering accuracy: 1.0000
Shower primary prediction accuracy: 0.0000
Track fragment clustering accuracy: 1.0000
Interaction grouping accuracy: 0.9000
Particle ID accuracy: 1.0000
Primary particle score accuracy: 0.8000
collect_dqdx
array([ 8383.44787459, 3075.15671474, 856.18962432, 10275.39597007,
995.6133342 , 1176.19442924, 5094.97991239, 1606.28092949,
2024.99650065, 1587.39534797, 2048.33477851, 7871.77062329,
6345.57253712, 2074.28885209, 1836.95889538, 3458.22282406,
6711.89453074, 4442.57608566, 5655.62416876, 2084.21310505,
15713.36874006, 3607.98016569, 4140.71172742, 2981.38425145,
3975.09149345, 3203.11402457, 4633.08825372, 1734.03140146,
3088.19267849, 2539.68349723, 2164.68916162, 2319.63803613,
1835.70240729, 718.86248216, 3659.63633432, 583.81144529,
4005.20169885, 7393.89333927, 5259.90855816, 2184.31200678,
1609.40608415, 4267.68260069, 3879.05482131, 5105.65571194,
6185.2047764 , 2061.0527234 , 1471.33420753, 1329.71173465,
5702.04481627, 3774.81705912, 6619.37141854, 3663.30712575,
1966.93197367, 2857.099469 , 4209.95435594, 2446.10420369,
3168.27279198, 1770.08554637, 2131.83559023, 3888.15439601,
1450.68970788, 3549.80349563, 5401.61544758, 1235.30559812,
4545.74656829, 1694.24682368, 1389.23629045, 5435.98994269,
3248.57105291, 3641.49973808, 4781.33406453, 1937.21630116,
3668.666896 , 2377.0897961 , 5269.51972236, 17243.41468401,
3612.28862708, 1283.01055095, 2295.38587576, 2151.81261546,
3672.65459286, 3226.7733467 , 4188.05570113, 4453.49283629,
5639.36500618, 4145.05137581, 4346.63280712, 1686.36702339,
11644.26216571, 1889.93101748, 1511.18817736, 6765.989624 ,
2818.90808417, 3852.3753254 , 4453.06565724, 3311.37873265,
2342.01964364, 4185.40556866, 1791.47230595, 1743.04679476,
4944.31346222, 11237.37441308, 1555.4987397 , 4080.65732143,
6638.5093132 , 2497.58095871, 2450.05506897, 1608.99026802,
2029.5626435 , 4258.23008719, 1635.99770029, 1473.70338119,
7415.35093279, 4242.0805526 , 1660.48735056, 1260.2668119 ,
1791.7354059 , 1753.12359994, 1931.13545477, 3858.60545567,
3266.70896608, 2716.79877939, 1330.92443848, 4916.78392715,
4648.18248606, 1604.38325131, 2205.23421692, 1997.71039713,
4825.46291631, 1463.48733788, 5612.17482923, 7909.49298996,
1278.4333537 , 2794.14678953, 864.87376771, 7823.02980173,
4710.22722834, 2802.62751409, 1642.87452029, 2451.03696642,
2158.67235926, 4567.52157391, 1908.14601307, 1222.02570602,
1349.19917312, 2261.88207952, 4567.61153747, 1629.47362839,
8685.72533064, 5874.09140703, 5191.60880059, 1592.21128529,
1113.48277196, 3108.80085299, 1606.12466419, 2898.62098508,
4306.86633509, 4433.66329358, 1303.48885931, 2065.59814292,
4443.62302853, 1989.77595656, 2768.49099097, 7455.54172597,
3542.01472489, 2393.80506845, 2191.05456677, 4379.83122454,
1344.33112751, 1940.17378688, 2347.78723609, 4588.90142353,
2565.57427164, 1625.41132676, 5078.9218029 , 4680.11439817,
4145.24923224, 5776.27695821, 1585.2462126 , 2839.11812281,
1462.0532709 , 3156.73323287, 2813.7863439 , 2724.78354291,
3591.54804823, 2079.11549178, 2481.55748291, 6557.37032988,
2664.00007719, 8274.65888803, 4066.90042351, 2021.37760225,
2626.39918086, 1706.91480243, 3266.29862745, 6556.13710787,
1612.54121498, 1917.43368685, 3339.60850803, 6588.74190669,
1799.6944859 , 2897.46137901, 2124.82437166, 1860.61781809,
2397.91343631, 3105.56355578, 3786.0666232 , 2374.60806139,
5266.3662759 , 6267.56842647, 7004.7296959 , 980.07165254,
5866.29261028, 8216.50239852, 2186.55955747, 2037.67489369,
7229.51468572, 1794.39217209, 11365.73869412, 3881.46629047,
5523.23893662, 5056.58824364, 1973.44170723, 2830.31733061,
1681.02516479, 6690.00966216, 2277.16182534, 2123.06629568,
2232.60620872, 2740.20779195, 4286.36090996, 1906.32678771,
1934.76181214, 2081.66136553, 4825.26648769, 3426.94249517,
3139.10820024, 1923.43318507, 1722.74570783, 4634.504832 ])
import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
'figure.figsize':(15, 10),
})
seaborn.set_context('talk')
plt.hist(collect_dqdx, range=[0, 10000], bins=50)
plt.xlabel("dQ/dx")
plt.ylabel("Predicted primary shower fragments")
Text(0, 0.5, 'Predicted primary shower fragments')
